Thermal stimulus task fMRI in the cervical spinal cord at 7 Tesla

Abstract Although functional magnetic resonance imaging (fMRI) is widely applied in the brain, fMRI of the spinal cord is more technically demanding. Proximity to the vertebral column and lungs results in strong spatial inhomogeneity and temporal fluctuations in B0. Increasing field strength enables higher spatial resolution and improved sensitivity to blood oxygenation level‐dependent (BOLD) signal, but amplifies the effects of B0 inhomogeneity. In this work, we present the first task fMRI in the spinal cord at 7 T. Further, we compare the performance of single‐shot and multi‐shot 2D echo‐planar imaging (EPI) protocols, which differ in sensitivity to spatial and temporal B0 inhomogeneity. The cervical spinal cords of 11 healthy volunteers were scanned at 7 T using single‐shot 2D EPI at 0.75 mm in‐plane resolution and multi‐shot 2D EPI at 0.75 and 0.6 mm in‐plane resolutions. All protocols used 3 mm slice thickness. For each protocol, the BOLD response to 13 10‐s noxious thermal stimuli applied to the right thumb was acquired in a 10‐min fMRI run. Image quality, temporal signal to noise ratio (SNR), and BOLD activation (percent signal change and z‐stat) at both individual‐ and group‐level were evaluated between the protocols. Temporal SNR was highest in single‐shot and multi‐shot 0.75 mm protocols. In group‐level analyses, activation clusters appeared in all protocols in the ipsilateral dorsal quadrant at the expected C6 neurological level. In individual‐level analyses, activation clusters at the expected level were detected in some, but not all subjects and protocols. Single‐shot 0.75 mm generally produced the highest mean z‐statistic, while multi‐shot 0.60 mm produced the best‐localized activation clusters and the least geometric distortion. Larger than expected within‐subject segmental variation of BOLD activation along the cord was observed. Group‐level sensory task fMRI of the cervical spinal cord is feasible at 7 T with single‐shot or multi‐shot EPI. The best choice of protocol will likely depend on the relative importance of sensitivity to activation versus spatial localization of activation for a given experiment. Practitioner Points First stimulus task fMRI results in the spinal cord at 7 T. Single‐shot 0.75 mm 2D EPI produced the highest mean z‐statistic. Multi‐shot 0.60 mm 2D EPI provided the best‐localized activation and least distortion.

generally produced the highest mean z-statistic, while multi-shot 0.60 mm produced the best-localized activation clusters and the least geometric distortion.Larger than expected within-subject segmental variation of BOLD activation along the cord was observed.Group-level sensory task fMRI of the cervical spinal cord is feasible at 7 T with single-shot or multi-shot EPI.The best choice of protocol will likely depend on the relative importance of sensitivity to activation versus spatial localization of activation for a given experiment.

Practitioner Points
1. First stimulus task fMRI results in the spinal cord at 7 T.
3. Multi-shot 0.60 mm 2D EPI provided the best-localized activation and least distortion.

| INTRODUCTION
The spinal cord houses many scientifically and clinically important neural circuits, such as central pattern generators related to locomotion (Minassian et al., 2017), sensorimotor reflex arcs, and circuits involved in the modulation of nociception (Bingel & Tracey, 2008).
Although functional magnetic resonance imaging (fMRI) is a wellestablished and widely applied imaging modality to probe neural activity in the brain, fMRI of the spinal cord is more technically demanding.
The adult human spinal cord is approximately 15 mm in diameter at its widest point, with gray matter structures only approximately 1-4 mm wide.High spatial resolution is therefore required to image the blood oxygenation level-dependent (BOLD) signal in the spinal cord.
In addition to its small size, several additional factors complicate fMRI of the spinal cord.The spinal cord is directly adjacent to the vertebral column, which is composed of alternating vertebral bones and soft tissues that differ in magnetic susceptibility, introducing distortions in the B 0 field.These spatially periodic B 0 inhomogeneities in the spinal cord cannot be adequately canceled by a counter-field generated by the scanner's shim coils, leaving residual focal B 0 inhomogeneities near the intervertebral junctions.The spinal cord is also near the lungs, whose volume and concentration of paramagnetic O 2 varies over time due to respiration, also introducing B 0 distortions.The resulting temporal fluctuations in B 0 in the spinal cord (Raj et al., 2001) produce frequency variation that can exceed 100 Hz at C7 (Vannesjo et al., 2017).These spatially and temporally varying B 0 inhomogeneities produce geometric distortion and signal loss, as well as apparent motion between volumes in the time series.Finally, pulsatile rostral-caudal flow of cerebrospinal fluid (CSF) over the cardiac cycle, and intermittent swallowing and motion, can also contribute nuisance variations to fMRI data that are more extreme than in the cerebrum (Cohen-Adad, 2017; Krüger & Glover, 2001).
Increasing the static magnetic field (B 0 ) improves signal to noise ratio (SNR), which translates to higher achievable spatial resolution.
Increased B 0 also increases the BOLD signal acquired by gradientecho (GRE) EPI, leading to better sensitivity and spatial specificity (Gati et al., 1997;Yacoub et al., 2001).These advantages motivate the use of 7 T fMRI to investigate circuits in the spinal cord.However, several difficulties are magnified and become significant at 7 T or higher field strengths.As field strength increases, susceptibilitymediated effects are magnified.Hence, although SNR and BOLD signal are enhanced, the geometric distortions and signal loss associated with B 0 inhomogeneity arising from the lungs and vertebral column are also intensified (Stroman et al., 2014;Triantafyllou et al., 2005).
Most 7 T scanners also lack a body RF transmit coil, and standard head RF coils do not cover the cervical spinal cord, which necessitates specialized RF hardware for the spinal cord that remains only accessible at a handful of research centers (Barry et al., 2018).
These challenges have hampered the development of spinal cord fMRI at 7 T, despite the potential advantages of higher field strength.
To date there are only three studies of spinal cord fMRI at 7 T (Barry et al., 2014(Barry et al., , 2016;;Conrad et al., 2018) all of which are resting-state fMRI studies investigating connectivity between gray matter regions in the absence of external stimuli or motor tasks.A significant motivation for the development of spinal cord fMRI methods, however, is to study pain processing in healthy and pathological states (Brooks & Tracey, 2005;Eippert et al., 2009;Geuter & Büchel, 2013;Staud et al., 2021;Tinnermann et al., 2017;Weber et al., 2016a).This involves measuring the BOLD response to external painful stimuli.
There is thus a need for investigating the potential of 7 T for taskbased spinal cord fMRI.
One important step toward fully realizing the potential benefits of 7 T is to establish suitable acquisition sequences and protocols that yield high BOLD sensitivity while being sufficiently robust against the challenges posed by static and dynamic B 0 inhomogeneities.At lower field strengths, various sequences have been investigated for spinal cord fMRI.Single-shot echo-planar imaging (EPI) is currently the dominant method for fMRI in the human spinal cord, as in the brain.The initial motor task fMRI studies in the cervical spinal cord, however, used a variety of pulse sequences, including fast low-angle shot (FLASH) sequences at 1.5 T (Yoshizawa et al., 1996) and 3 T (Stroman et al., 1999), multi-shot gradient echo (GRE) EPI at 1.5 T (Backes et al., 2001), spin-echo EPI at 1.5 T (Stroman & Ryner, 2001), and single-shot fast spin echo at 1.5 T (Stroman, 2002;Stroman et al., 2002a).Spin-echo-based pulse sequences continued to be used for fMRI in the spinal cord at 1.5 T based on signal enhancement by extravascular water protons (SEEP) contrast for many years (Lawrence et al., 2008;Stroman et al., 2002b).Work on BOLD fMRI in the spinal cord by multiple groups over the next several years focused mainly on 2D single-shot (Brooks et al., 2008;Eippert et al., 2016;Giulietti et al., 2008;Govers et al., 2007;Piché et al., 2009;Weber et al., 2016bWeber et al., , 2020) ) and 3D multi-shot (Barry et al., 2020) GRE EPI at 1.5 T and 3 T.
Most relevant to the present work, others have previously implemented 2D single-shot EPI in combination with noxious thermal stimulation to study cognitive manipulation of pain (Eippert et al., 2009;Sprenger et al., 2012) the effects of noxious thermal stimulation on functional connectivity in the cervical spinal cord (Weber et al., 2017), and several have developed and implemented methods for simultaneous brain and spinal cord MRI (Chu et al., 2023;Finsterbusch et al., 2013;Islam et al., 2018) and utilized them in studies of, for example, the fundamental organization of resting-state networks spanning the brain and the cervical spinal cord (Vahdat et al., 2015(Vahdat et al., , 2020) ) and the interaction between brain and spinal cord regions on sensitivity to pain (Sprenger et al., 2015;Tinnermann et al., 2017).
Single-shot EPI sequences are heavily affected by geometric distortion and signal loss caused by static B 0 inhomogeneity.Multi-shot EPI sequences have strong potential to reduce these artifacts by shortening the echo train length and echo time.In the previous 7 T resting-state studies, highly segmented multi-shot 3D EPI with very short echo train lengths was used (Barry et al., 2011(Barry et al., , 2014(Barry et al., , 2016)), yielding images largely uncorrupted by signal loss and distortion near intervertebral junctions.However, the transition from single-shot to multi-shot EPI entails a penalty in volume acquisition time (VAT), which can reduce statistical power to detect to BOLD activation by reducing the number of observation time points for a given scan duration.Multi-shot EPI is also more vulnerable to dynamic B 0 field fluctuations, resulting in variable ghosting between volumes in the time series, which can reduce the temporal SNR (Vannesjo et al., 2019).In single-shot EPI, on the other hand, time-variable B 0 fields (at the respiratory time-scale) largely translate into apparent motion in the images, which can be addressed by retrospective motion correction.
Considering the relative advantages and disadvantages of single-shot versus multi-shot EPI sequences, it is not apparent which is more suitable for fMRI of the spinal cord at ultra-high field.
In this work, we present the first stimulus task fMRI in the cervical spinal cord at 7 T and compare the performance of one single-shot 2D EPI and two multi-shot 2D EPI protocols in terms of sensitivity and specificity to stimulus-related activation.

| Image acquisition
The cervical spinal cords (C4-C7 vertebral levels, approximately corresponding to spinal cord segments C5-C8) of 11 healthy volunteers (ages 22-39, 5 males, 6 females) were scanned using a 7-T wholebody MRI system (Magnetom 7 T AS, Siemens Healthineers, Erlangen, Germany) equipped with a G max = 70 mT/m and maximum slew rate = 200 mT/m/ms gradient coil and a single-channel four-element transmit, 22-channel receive array brainstem and cervical spine radiofrequency (RF) coil (Zhang et al., 2017) One single-shot 2D EPI protocol at 0.75 mm in-plane resolution and two multi-shot 2D EPI protocols at 0.75 and 0.6 mm in-plane resolution were used (Figure 1).These three protocols will be abbreviated SS75, MS75, and MS60, respectively.Each protocol was optimized independently, rather than attempting to match parameters between protocols.Research pulse sequences were used for both the single-shot and multi-shot acquisitions.(At the time of this study, a multi-shot 3D EPI pulse sequence analogous to what was used in ref- erences [Barry et al., 2014[Barry et al., , 2016;;Conrad et al., 2018] was not available on our scanner platform.)In all protocols, the field of view (FOV) was intentionally centered toward the back of the neck to include posterior neck muscles, which, due to their greater distance from the trachea and esophagus, serve as spatial landmarks that are relatively robust to swallowing to aid subsequent motion correction.
In the single-shot EPI, a single saturation slab was placed anterior to the spinal cord for outer volume suppression, thereby reducing the required echo-train length.For in-plane acceleration using generalized autocalibrating partially parallel acquisition (GRAPPA) (Griswold et al., 2002), 52 lines of autocalibration signal reference data were acquired using GRE (Talagala et al., 2016), which has been shown to be optimal for single-shot EPI in the spinal cord at 7 T (Seifert & Xu, 2022).Partial Fourier (factor 6/8) acquisition was judged to be necessary in single-shot EPI to further reduce echo time (TE) to avoid excessive signal loss due to through-slice dephasing.Image reconstruction of the single-shot EPI was performed online with the standard reconstruction pipeline of the manufacturer, including partial Fourier reconstruction by zero-filling and sine-window filtering, and Nyquist ghost correction based on a three-line navigator.
The multi-shot protocols were acquired without a spatial saturation slab due to specific absorption rate (SAR) limitations given the shorter per-shot repetition time (TR).Extending the VAT to accommodate spatial saturation was judged to be an unfavorable tradeoff.
Therefore, a larger FOV had to be acquired in the multi-shot protocols, than in the single-shot protocol, to encompass signal from the entire neck.Note that phase-encode FOVs cannot be defined to be larger than readout-encode FOV, and the readout-encode FOV can be enlarged with no SNR or timing penalties.For this reason, the FOV in the multi-shot protocol was set larger in the readout (right-left) dimension than in the single-shot protocol to achieve the intended increase in the phase-encode FOV.Phase-encoding in the left-right dimension is undesirable because the patterns of spatial distortion in the spinal cord would be asymmetric across the left-right dimension of the spinal cord.Unlike the single-shot sequence, the custom multishot sequence did not allow for GRAPPA acceleration via the console, therefore the FOV in the phase-encode dimension was set to half of the readout FOV in the multi-shot protocols to achieve the intended phase-encode undersampling (Figure 1).Image reconstruction of the multi-shot acquisitions was performed offline, with an iterative conjugate gradient SENSE algorithm (Pruessmann et al., 2001).The coil sensitivities were estimated based on a low-resolution GRE acquisition.
For SNR-optimal coil combination, the noise correlation between coils was calculated based on data acquired in vivo without RF excitation.
A three-line navigator was acquired for each shot, which was used for per-shot frequency offset demodulation to reduce ghosting caused by temporal B 0 field fluctuations (Vannesjo et al., 2019).The frequency offset per shot was estimated based on a magnitude-weighted coil combination of the phase at the center of k-space in the second navigator echo.For Nyquist ghost correction, three-line phase reference data (which characterize the gradient waveforms) were also acquired in a phantom immediately following each session, using the same multi-shot protocols as performed in vivo.In pilot studies, this approach was observed to be more effective in reducing Nyquist ghosting than using the three-line navigators from the in vivo data.
To compare expected advantages and disadvantages between the protocols, Table 1 lists performance considerations relevant to BOLD fMRI, and for each of these gives the governing protocol parameters, and the expected relative ranking of the three candidate protocols.
The relative importance of the different performance considerations, and the holistic ranking of protocol performance in cervical spinal cord BOLD fMRI at 7 T, is not obvious; the goal of this study is to determine this experimentally.

| Functional MRI task design
For each of the three imaging protocols, a 10-min fMRI acquisition was performed with noxious thermal stimulation applied to the lateral surface of the base of the right thumb with a 30 mm Â 30 mm flat thermode at a pre-calibrated intensity of 3/10 (range 42-51 C, average 46.5 C) using an fMRI-compatible thermal stimulator (TSA-II, Medoc, Israel).Pain scale anchors were 0 = no pain and 10 = worst pain imaginable.Stimulus blocks were 10 s in duration, calculated by including ramp-up and plateau, but excluding ramp-down to baseline.
Ramp speed was set to 10 C/s.Stimulus blocks were separated by pseudo-randomized inter-stimulus intervals ranging from 25 to 45 s, plus initial and final stimulus-free blocks.This stimulus was expected to produce sensory activation primarily in the ipsilateral dorsal horn at the neurological C6 (approximately vertebral C4-C5) level (Keegan & Garrett, 1948;Lee et al., 2008).Subjects were coached to remain still and maintain steady respiration during thermal stimuli.Pulse-oximeter and respiratory traces were acquired (TSD200, PPG100C, TSD201, and RSP100C, Biopac, Goleta, CA, USA).

| Data processing
Spinal cord toolbox (SCT) v5.4.0 (De Leener, Lévy, et al., 2017) was used for image pre-processing.Anatomical MP2RAGE images were segmented using sct_deepseg_sc, and vertebral levels were labeled using sct_label_vertebrae and/or sct_label_utils.Registration to the PAM50 template (De Leener, Fonov, et al., 2017) was performed using sct_register_to_template.Each 4-D fMRI dataset was motioncorrected using sct_fmri_moco with an inclusive cylindrical mask that contains the posterior neck muscle, the spinal canal and spinal cord were manually masked in FSLeyes (Jenkinson et al., 2012), and the fMRI dataset was registered to the PAM50 template using sct_regis-ter_multimodal, initialized using the transformation calculated between the anatomical image and PAM50.Each fMRI dataset was straightened, spatially smoothed using sct_smooth_spinalcord with an anisotropic 2 mm Â 2 mm Â 6 mm Gaussian kernel oriented along the spinal cord centerline, and unstraightened, to minimize blurring of CSF signal into the cord.
First-level general linear model (GLM) analysis was performed in individual-subject space using FSL FEAT (Woolrich et al., 2001), incorporating the default autocorrelation correction using FSL FILM (Woolrich et al., 2001) and a 35-term physiological noise model (PNM) (Brooks et al., 2008;Kong et al., 2012).The PNM was generated separately for each slice, incorporating slice timing correction, using FSL PNM, and included 8 cardiac phase, 8 respiratory phase, 16 cardiac-respiratory interaction, CSF signal derived from the spinal canal mask (canal mask minus cord mask), and two in-plane translational motion correction terms from sct_fmri_moco.
Group-level GLM analysis was performed in FSL FEAT using a mixed-effects model (FLAME1).After warping the results of first-level analyses into PAM50 template space, group-level activation maps were calculated within a search region constructed by eroding the PAM50 cord mask by one voxel.First, we performed thresholding at z > 2.7 followed by cluster-level correction at a cluster threshold of p = .05.This cluster-level correction method was tested because of its widespread use and familiarity in the field of brain fMRI, but has not yet been rigorously validated in the spinal cord, and so these results should be interpreted with caution.Next, because cluster-level correction did not yield any supra-threshold activation for any of the three protocols, we applied a low, uncorrected threshold of z > 1.64, which corresponds to a voxel p-value threshold of 0.05.This uncorrected threshold of z > 1.64 was subsequently used for tabulation of activated voxel counts and visualization in support of the goal of comparing the three acquisition protocols.

| Data analysis
Temporal signal to noise ratio (tSNR) was calculated voxelwise as the ratio of temporal mean to standard deviation, then averaged within neurological level C6 in the PAM50 cord, gray matter, and white matter masks after motion correction, but before smoothing and application of the PNM.The extent of this neurological C6 mask was defined generously as the region where the probabilistic mask of C6 contained within SCT (Cadotte et al., 2015) is greater than zero.Several additional region of interest (ROI) masks were defined for further analysis ).These masks were also warped to each individual subject space.Although studies in rats have shown that some contralateral ventral horn activation does occur in response to noxious stimuli (Coghill et al., 1991(Coghill et al., , 1993)), this was judged to be the most suitable control region in light of the stronger resting-state connectivity between ipsilateral dorsal and ventral horns, and between ipsilateral and contralateral dorsal horns, in humans (Barry et al., 2014(Barry et al., , 2016;;Kong et al., 2014).
For the group-level and each individual-subject result, the mean and maximum z-statistic across voxels with z > 1.64 within

| RESULTS
Individual subject-level mean and single-timeframe images and tSNR maps of datasets representing the best-case and worst-case image quality in the study are shown in Figure 3.These best-case and worst-case datasets were chosen by first selecting the subsets of cases with the highest and lowest tSNR, and from those subsets, choosing the dataset with the lowest and highest artifact burden by visual inspection, respectively.In axial slices near intervertebral junctions, where residual B 0 inhomogeneity is large, signal loss due to through-slice dephasing and geometric distortion are evident.Geometric distortion was less severe in multi-shot acquisitions, which had higher effective phase encoding bandwidth.The single-shot EPI had a blurrier appearance than the multi-shot acquisitions, while multi-shot 0.60 mm produced the sharpest appearing images.A spatial gradient in tSNR was also evident, with lower tSNR at more inferior levels.This is consistent with a combination of reduced B 1 + efficiency, increased static B 0 inhomogeneity, and larger time-varying B 0 fluctuations due to the proximity to the lungs.
Summary statistics of tSNR are given in Table 2. Temporal SNR was generally highest in single-shot 0.75 mm and multi-shot 0.75 mm, with lower variability between subjects in multi-shot 0.75 mm, while tSNR was significantly lower in multi-shot 0.6 mm, yet with lowest between-subject variability.Temporal SNR was slightly greater in gray matter than in white matter in all three protocols, driven predominantly by the higher signal in gray matter.
Using cluster-level correction for multiple comparisons within the eroded PAM50 cord mask, no significantly activated voxels were identified in the spinal cord on group level analyses for any of the three  3. Voxels within the true-positive and control ROI masks with z > 1.64 are considered "activated" for the purpose of comparisons between protocols, although as stated previously, these tabulated voxels do not survive multiple-comparison correction using cluster-level correction.Mean z-statistics in activated true-positive voxels are relatively consistent in the range of 1.83-2.04across the three protocols, with single-shot 0.75 mm yielding the highest mean z-statistic contrary to initial expectation (Table 1).The maximum truepositive z-statistic is highest in single-shot 0.75 mm at 2.85, compared to 2.07 in multi-shot 0.75 mm and 2.38 in multi-shot 0.6 mm.All three protocols yielded a ratio of true-positive to control max z-statistic between 1.27 and 1.37.
The mean BOLD PSC in activated voxels was higher in multi-shot datasets than in single-shot.Multi-shot 0.75 mm yielded the highest mean BOLD PSC, at 0.43%, and the BOLD PSC in the voxel with the highest z-statistic was 0.50%.
F I G U R E 2 Masks used for analysis of task fMRI activation, displayed in PAM50 template space.The PAM50 cord mask, eroded by one voxel (blue), was used as the search region in group-level analyses, and to mask all parameter maps for visualization.Masks representing regions of interest in the C6 spinal level for expected true positive activation (green) and control (red) are also displayed.
To estimate the number of required subjects to achieve a suprathreshold cluster using stringent cluster-level correction methods, we performed additional group-level analyses with cluster-level correction as described in the methods section, on datasets containing replicates of the set of 11 individual subject-level datasets.Based on these analyses, we estimate that using thresholding at z > 2.7 followed by cluster-level correction at a cluster threshold of p = .05(which might be considered "common" methods in the brain), the required number of subjects is estimated to be between 23 and 33 in single-shot 0.75 mm, between 177 and 187 in multi-shot 0.75 mm, and between 34 and 44 in multi-shot 0.60 mm.

Details of z-statistics and BOLD PSC across all individual-level
GLM analyses are given in Table 4.As in group-level analyses, voxels within the true-positive and control ROI masks (warped to each individual subject space) with z > 1.64 are considered activated.Mean zstatistics in activated true-positive voxels are relatively consistent in the range of 1.96-2.18across the three protocols.Single-shot 0.75 mm yields the highest mean z-statistic, but the standard deviation of each subject's mean z-statistic in activated true-positive voxels is double that of the multi-shot protocols.Single-shot 0.75 mm also yields the highest mean maximum z-statistic at 2.68, compared to 2.27 for multi-shot 0.75 mm and 2.45 for multi-shot 0.6 mm.Singleshot 0.75 mm and multi-shot 0.6 mm have true-positive-to-control ratios substantially greater than unity, while the ratio in multi-shot 0.75 mm is closer to unity.
Similar to group-level analyses, the mean BOLD PSC in activated voxels was higher in multi-shot protocols.Multi-shot 0.6 mm yielded the highest mean BOLD PSC, at 1.78 ± 0.89%, followed by multi-shot 0.75 mm at 1.29 ± 0.76%, and single-shot 0.75 mm at 0.67 ± 0.31%.
The same relationship among the protocols is evident in the BOLD PSC in the voxel with the highest z-statistic in each dataset.T A B L E 2 Summary statistics of temporal signal to noise ratio (tSNR) measurements in data from three protocols.Note: Tabulated means and standard deviations are calculated as mean and standard deviation of the set of individual subject-level means within PAM50 masks.NB: standard deviations across subjects of mean tSNR in the cord, GM, and WM do indeed all round to 4.2; differences appear at the third decimal place.
tSNR was sharply reduced and geometric distortion and signal loss were markedly more severe than in slices further from a junction.The gradual reduction in tSNR in lower slices was partially due to static and dynamic B 0 inhomogeneity from the lungs, but reduced B 1 + efficiency in the lower cervical spine was also a significant contributor (Zhang et al., 2017); signal intensity in mean and single images was lower adjacent to vertebral level C7, compared to C4-C6, even in the best-quality dataset.
F I G U R E 4 Maps of z-statistics and BOLD percent signal changes in grouplevel GLM analyses.z-Statistic maps are displayed using a low, uncorrected threshold of z > 1.64, and as such, should be interpreted with caution.Maps of BOLD percent signal change are thresholded at 0.125% for ease of visualization.Maps are masked to the spinal cord, encompassing all neurological levels, and the slice locations in all three planes are centered on the voxel inside the spinal cord mask with the highest zstatistic (indicated with a green crosshair).
T A B L E 3 Summary statistics of group-level GLM results.
Group-level activation results

SS75 MS75 MS60
Mean z-stat in activated voxels in true positive region 2.04 Note: An uncorrected threshold of z > 1.64 was used to tabulate activated voxel counts for the purpose of comparing acquisition protocols; therefore, these results should be interpreted with caution.
In general, the single-shot 0.75 mm protocol provided the greatest sensitivity (albeit highest between-subject variability in z-statistics) to the expected activation, but also appeared to have the lowest effective in-plane resolution and largest geometric distortion.This was in line with expectations (Table 1), as the lower effective phaseencode bandwidth makes it more susceptible to T 2 * blurring and geometric distortion due to static B 0 inhomogeneity, on top of blurring due to partial Fourier acquisition.Multi-shot 0.6 mm produced the most well-localized activation clusters, and the sharpest images, while being qualitatively less affected by geometric distortion caused by static B 0 inhomogeneity than single-shot 0.75 mm.However, sensitivity to activation was reduced due to its longer VAT, which reduced the number of time points in an experiment of fixed duration, and the lower tSNR, both of which reduce the statistical power.Multi-shot 0.75 mm generally performed poorly; it neither showed the sensitivity to activation of the single-shot 0.75 mm protocol nor the spatial specificity of the multi-shot 0.60 mm protocol, as if sharing the drawbacks of each of the other protocols, without enjoying the benefits of either protocol to a similar degree.
Decades of results in the brain have consistently shown that longer TEs (i.e., close to the T 2 * of gray matter) are associated with greater sensitivity to activation.In the spinal cord at 7 T, however, the detrimental effects of B 0 inhomogeneity severely impair data quality at longer TEs.A balance must be struck between image quality, which favors shorter TE, and sensitivity to BOLD signal, which favors TE close to T 2 * of gray matter, which is 29.3 ± 4.5 ms at 7 T (Massire et al., 2016).The optimal TE for a given spinal cord fMRI pulse sequence will therefore likely be shorter than the T 2 * of gray matter.
In group-level GLM analyses, single-shot 0.75 mm and multi-shot 0.6 mm appeared to yield similarly sized activation clusters, and the maximum z-statistic in single-shot 0.75 mm was substantially greater than the maximum z-statistic in multi-shot 0.6 mm (see Figure 4).In individual-level results, however, multi-shot 0.6 mm generally yielded more sharply defined activation clusters in the ipsilateral dorsal quadrant (see Figure S1), and the difference between the maximum zstatistic in single-shot 0.75 mm and the maximum z-statistic in multishot 0.6 mm was much smaller than in the group level analysis.A possible explanation that unifies these observations is that multi-shot 0.6 mm was better able to spatially resolve activation in individual subjects, resulting in small clusters; however, when multiple subjects' data was combined, the inherent physiological variability in the location of activation among different subjects (Cadotte et al., 2015), as well as potential limitations in group-level registration of distorted spinal cord EPI images, resulted in these sharply defined individual-level activation clusters spatially coinciding less fully in multi-shot 0.6 mm than did the broader activation clusters in single-shot 0.75 mm.These broader activation clusters in single-shot EPI may have arisen due to its greater statistical power, having more temporal frames due to a shorter VAT and higher tSNR due to the lower effective resolution.
As stated previously, the noxious thermal stimulus was expected to produce sensory activation in the ipsilateral dorsal horn at the neurological C6 (approximately vertebral C4-C5) level based on the fact that the majority of first-order neurons in the lateral spinothalamic tract synapse in the ipsilateral dorsal horn (Keegan & Garrett, 1948;Lee et al., 2008).Indeed, virtually all previous spinal cord fMRI studies of noxious thermal stimulation at lower field strengths or in animals T A B L E 4 Summary statistics of individual subject GLM analyses.
Individual-level activation results

SS75 MS75 MS60
Mean z-stat in activated voxels in true positive region 2.18 ± 0.37 1.99 ± 0.17 have observed the strongest activation in the ipsilateral dorsal horn.
Based on TE, BOLD PSC would be expected to be highest in multi-shot 0.60 mm, followed by single-shot 0.75 mm and then multishot 0.75 mm.In individual-level results, multi-shot 0.60 mm did indeed yield the highest BOLD PSC.However, of the remaining two protocols, multi-shot 0.75 mm yielded higher BOLD PSC than singleshot 0.75 mm.One potential explanation for this is that the effect of spatial resolution outweighed the effect of the magnitude of the underlying BOLD activation in these individual-level measurements; the multi-shot 0.75 mm acquisition is expected to better resolve BOLD activation than the single-shot 0.75 mm acquisition.An alternative explanation specifically concerning the difference in the mean BOLD PSC in activated voxels between the multi-shot 0.75 mm and single-shot 0.75 mm acquisitions is that the single-shot 0.75 mm acquisition, which contained 400 temporal frames, was statistically better powered to detect low-BOLD PSC activation than the multishot 0.75 mm acquisition, which contained only 202 temporal frames and may have been additionally hampered by greater physiologic noise contributions (e.g., ghosting).In other words, a voxel with a low BOLD PSC is more likely to be identified as significantly activated using a protocol with a larger number of temporal frames and less susceptibility to ghosting.This would result in a greater number of low-BOLD PSC voxels being included in the set of significantly activated voxels in the single-shot 0.75 mm acquisition, decreasing the mean BOLD PSC calculated for that set of voxels.However, when the mean BOLD PSC is calculated for each protocol in the set of voxels having z-statistics in the top 10th percentile within each individual subject, BOLD PSC remains lower in the single-shot 0.75 mm acquisition than in the multi-shot 0.75 mm acquisition.In the same percentile-based set of voxels in group-level results, single-shot 0.75 mm likewise has the lowest mean BOLD PSC.In light of these results, the first potential explanation based on spatial resolution may be more likely.
The superior-inferior location of peak activation was similar across the three protocols in group-level analyses, and appeared at the expected location in the center of neurological C6 (determined from the probabilistic spinal levels from the PAM50 atlas), based on the anatomical location of the thermal stimulus.The superior-inferior location of the peak activation was expected to vary across subjects due to natural physiological variation (Cadotte et al., 2015), but was not expected to vary among the three protocols within each subject.
Contrary to this expectation, we did observe substantial variation in the superior-inferior location of peak activation of up to 1½ vertebral levels across the three protocols within individual subjects (see Figure S1).The cause for this variation is unclear without conducting a reproducibility study by repeating all three protocols on the same participants.If the true underlying stimulus-related activation occurs across a large spatial extent in the superior-inferior dimension, then random noise as well as variations in physiological noise may shift the center of the detected activation between acquisitions in an unsystematic way, as observed.In that case, increasing the sensitivity to detect BOLD responses, for example, by acquiring for longer, should yield extended activation clusters.The true underlying activation is however unknown, and the data in this work cannot convincingly support any conclusion.From resting-state studies, we would expect networks of correlated activity stretching over approximately one vertebral level (Kong et al., 2014).Task-activated networks, and patterns of nociceptive activation in particular may, however, have a broader spatial extent to resting-state networks due to the fact that primary afferent projections may ascend or descend in Lissauer's tract before synapsing (Coggeshall et al., 1981;Light & Perl, 1979;Sugiura et al., 1986).A large spatial extent of the activated networks in the superior-inferior dimension may even be a prerequisite to detect activation clusters in group-level analyses with relatively few subjects.If the activation were narrowly localized within one subject, with a variable center along the vertebral levels between subjects, the statistical power to detect activation at the group level would be much diluted.
Furthermore, it has been shown in the cerebral cortex that significant variation in activation maps exists across sessions, even in a single subject (McGonigle et al., 2000), and when evaluating activation based on thresholded single-subject activation maps, even consistent activation can appear highly variable when the significance of this activation is near the defined threshold for visualization (Smith et al., 2005).Further studies would be needed to investigate the true extent of the activation, and the reproducibility of individual-level activation clusters.
Several groups have developed methods to address distortion and through-slice dephasing in single-shot GRE EPI at 3 T through slicewise dynamic shimming methods using either slice-specific gradient offsets (Islam et al., 2018) or z-gradient refocusing blips (Finsterbusch et al., 2012;Kaptan et al., 2022).These methods, applied to singleshot EPI, may be promising alternatives to multi-shot 2D or 3D EPI at 7 T but, at the time of this study, were not sufficiently refined to implement routinely at 7 T. Also 3D multi-shot EPI sequences suitable for implementation in fMRI experiments were not readily available at the time of this study, despite encouraging results demonstrating their robustness to B 0 inhomogeneity in resting-state fMRI studies at 7 T (Barry et al., 2014(Barry et al., , 2016)).Such sequences remain a highly promising option worthy of future investigation.This study also used different reconstruction methods for single-shot and multi-shot EPI: on-scanner GRAPPA and offline SENSE, respectively.An ideal solution would be to reconstruct single-shot EPI data offline using the same SENSE method as the multi-shot EPI data; however, due to data size and transfer rate limitations, we were not able to save raw data for singleshot fMRI for all subjects.
This study is limited by its small sample size of 11 healthy subjects.Because of this, no significantly activated voxels were identified in the spinal cord on group level analyses using cluster-level correction for multiple comparisons.To enable visual and quantitative comparisons between protocols, maps and tables used a low, uncorrected threshold of z > 1.64.Future studies in a larger number of subjects would be needed to achieve significance under appropriately stringent thresholding methods.
It should also be noted that the BOLD PSC maps displayed in subjects have previously been presented as conference abstracts (Seifert & Vannesjo, 2020, 2022).

1
Pulse sequence parameters and slice positioning for three fMRI protocols.Yellow boxes indicate image acquisition fields of view, and green boxes indicate frequency and shim adjustment volumes.T A B L E 1 Performance considerations relevant to BOLD fMRI, the protocol parameters governing performance regarding each image quality or sensitivity consideration, and the expected performance of the three candidate protocols relative to each other.The bolded protocol names are the protocols that are expected to perform best in each performance criterion.

(
Figure2).The PAM50 cord mask, eroded by one voxel, was used as the search region in group-level analyses and to mask parameter maps for visualization.A mask for expected true-positive activation was constructed by summing the PAM50 masks for the right dorsal horn, right fasciculus gracilis, right fasciculus cuneatus, and right lateral corticospinal tract (PAM50 atlas mask entries 1, 3, 5, and 35), binarizing, eroding by one voxel from the cord-CSF border only, and calculating the intersection of the result with the aforementioned PAM50 mask for spinal (i.e., neurological) level C6 thresholded at zero.White matter structures were included to enclose potential BOLD signal in veins draining the expected gray matter site of activation(Seifert & Vannesjo, 2020), but the spinal venous plexus on the exterior of the spinal cord, though expected to also exhibit BOLD activation(Cohen- Adad et al., 2010), was excluded due to the strong effects of CSF signal that cannot be fully accounted for by a PNM and the resulting strong potential for false-positive activation.A control mask, for comparison purposes, was calculated similarly using PAM50 masks for the left lateral corticospinal tract, left rubrospinal tract, left spinothalamic and spinoreticular tracts, and left ventral horn (PAM50 mask entries4, 8, 12, and 30 the true-positive and control ROI masks were tabulated.The ratio of the maximum z-statistic in the true-positive ROI mask to the maximum z-statistic in the control ROI mask in each subject was taken as a measure of the specificity of each protocol to truepositive activation.Mean BOLD percent signal change (PSC) was calculated in activated native space voxels in the true-positive region for each subject, and the BOLD PSC in the voxel in the truepositive region having the highest z-statistic in each subject was also tabulated.
protocols.Maps of z-statistics, thresholded at an uncorrected threshold of z > 1.64, and BOLD PSCs in group-level GLM analyses are shown in Figure 4. Clusters of activation under this uncorrected threshold appear in all three protocols in the ipsilateral (right) dorsal quadrant at the expected C6 neurological level, which is centered approximately at the C4-C5 intervertebral disc (see Figure 2 for superior-inferior extent).All images are displayed in radiological orientation convention.The ipsilateral (right) activation appears in the same cross-sectional location in single-shot 0.75 mm and multi-shot 0.6 mm; however, this cluster is weaker and located more peripherally in multi-shot 0.75 mm.In single-shot 0.75 mm, a cluster also appears in the contralateral (left) dorsal horn.Details of z-statistics and BOLD PSC in group-level GLM analyses are given in Table Image quality and tSNR are strongly impaired in spinal cord fMRI by static focal B 0 inhomogeneity near the intervertebral junctions and temporally variable B 0 inhomogeneity arising from the lungs.These two sources of B 0 inhomogeneity were clearly evident in the images and tSNR maps shown in Figure3.Near the intervertebral junctions, F I G U R E 3 Individual subject-level images and temporal signal to noise ratio (tSNR) maps spanning C4-C7 vertebral levels in datasets representing best and worst image quality in the study.Temporal SNR was calculated after motion correction, but before spatial smoothing or application of the physiological noise model.Locations of the three axial slices are indicated with yellow lines in the sagittal frames.The first and third displayed axial slices are superior and inferior to an intervertebral disc, and the second axial slice is located at an intervertebral disc.

Figure 4 ,
Figure 4, Figures S1 and S2 display the amplitude of the BOLD signal but do not provide information on the variability of the BOLD signal, and so they should be interpreted with appropriate caution, and in combination with the accompanying activation z-score maps.
An uncorrected threshold of z > 1.64 was used to tabulate activated voxel counts for the purpose of comparing acquisition protocols; therefore, these results should be interpreted with caution.In individual subjects, mean z-statistics and mean BOLD percent signal change (PSC) are calculated within voxels having z > 1.64, as well as BOLD PSC in the voxel having the highest z-statistic and in the set of voxels having zstatistics in the top 10th percentile.Maximum z-statistics in true positive and control regions are also calculated in individual subjects without a threshold on z-statistic.Means and standard deviations across these individual-subject mean and maximum values are then calculated and tabulated here.